This R notebook conducts various tests using the method outlined in Toward automated extraction and characterization of scaling regions in dynamical systems.

# Dependencies
library(ggplot2)
library(LaplacesDemon)
source("./generic.R")
source("./deshMethod.R")
states_visited <- c()
acfs <- c()
# Specific to autoregression 
a = 0.95
b = 2
exp_m <- 1 + 2*(a/(1-a))
steps = 10000

get_next_step <- function(prev) {
  w <- rnorm(n=1, mean=0, sd=1)
  next_val <- a * prev + b * w
  return(next_val)
}

get_init_state <- function( ) {
  init_sd <- b^2 / (1-a^2)
  init_state <- rnorm(n=1, mean=0, sd=init_sd)
  return(init_state)
}
# example method with one simulation using one set of p, q
# also plots histograms here. 
n <- 10 # chosen
p <- 3 # chosen
q <- 3 # chosen 
steps = 10000
a <- 0.9
r <- 0

states_visited <- asim(get_init_state(), get_next_step, steps)
acf_df <- c_est_series(states_visited)

# limiting the search to a finite number, since takes a long time computationally
stop_search <- which(acf_df$acf <= 0)[1]
stop_search <- 1.5 * stop_search

y <- get_log(acf_df$acf)[1:stop_search]
x <- c(0:(stop_search-1))
lim <- length(y)
results_ar1 <- desh_method(x,y, lim, n)
[1] 58
Time difference of 0.2918248 secs
new_weights <- get_desh_weights(results_ar1,  p, q, r)
[1] 1128
vars_ar1 <- get_slopes_sides(results_ar1, new_weights, TRUE)
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

vars_ar1
exp_m <- 1 + 2*(a/(1-a))
sprintf("Resulting RHS : %f, expected time: %f", vars_ar1$best_rhs, exp_m)
[1] "Resulting RHS : 20.932509, expected time: 19.000000"

This is the big test– does the method for many different a’s, p’s, and q’s.

# Generates simulations for each possible a we want to test, storing results in a dictionary. 
#poss_as <- c(-0.9, -0.7, -0.5, 0.5, 0.7, 0.9, 0.99)
poss_as <- c(0.5, 0.7, 0.9, 0.99)
steps = 10000
dict <- c()
for(a in poss_as) {
  states_visited <- asim(get_init_state(), get_next_step, steps)
  acf_df <- c_est_series(states_visited)
  dict[as.character(a)] <- acf_df
}
# Then we try the method on each simulation.
# Start by finding the autocorrelation functions. 
poss_ps <- c(1, 2, 3)
poss_qs <- c(1, 2, 3)
poss_rs <- c(0)

old_ms <- c()
new_ms <- c()
ps <- c()
qs <- c()
as <- c()
real_m <- c()
results_ar1s <- list()

for(a in poss_as) {
  acf <- dict[as.character(a)]
  acf_df <- data.frame(acf)
  colnames(acf_df) <- c("acf")
  stop_search <- which(acf_df$acf <= 0)[1]
  stop_search <- 2 * stop_search
      
  if(a < 0) {
    n <- 3
  } else if(stop_search > 500) {
    stop_search <- 500
  } else {
    n <- 10
  }
  
  if(stop_search < 5) {
    stop_search <- 10
  }
  y <- get_log(acf_df$acf)[1:stop_search]
  x <- c(0:(stop_search-1))
  lim <- length(y)
  results_ar1 <- desh_method(x,y, lim, n)
  results_ar1s <- append(results_ar1s, results_ar1)
}
[1] 20
Time difference of 0.0153501 secs
[1] 36
Time difference of 0.1190212 secs
[1] 160
Time difference of 6.823408 secs
[1] 500
Time difference of 8.277671 mins
# Then try the method on the given autocorrelation functions.
old_ms <- c()
new_ms <- c()
ps <- c()
qs <- c()
as <- c()
real_m <- c()
rs <- c()

for(a in poss_as) {
  print(a)
  start_time <- Sys.time()
  c <- 4
  M_upper_bound <- 6000
  triangle_df <- find_m(acf_df, c, M_upper_bound)
  m_suggest <- triangle_df[triangle_df$poss_Ms >= triangle_df$ct_ints,]
  m_suggest <- m_suggest$poss_Ms[1]
  m_suggest <- t_int(m_suggest, acf_df)
  
  acf <- dict[as.character(a)]
  acf_df <- data.frame(acf)
  colnames(acf_df) <- c("acf")
  stop_search <- which(acf_df$acf <= 0)[1]
  stop_search <- 1.5 * stop_search
      
  if(a < 0) {
    n <- 3
  } else if(stop_search > 500) {
    stop_search <- 500
  } #else {
    #n <- 10
  #}
  
  if(stop_search < 5 ) {
    stop_search <- 10
  }
  y <- get_log(acf_df$acf)[1:stop_search]
  x <- c(0:(stop_search-1))
  lim <- length(y)
  results_ar1 <- desh_method(x,y, lim, n)
  
  
  for(p in poss_ps) {
    for(q in poss_qs) {
       for(r in poss_rs) {
        exp_m <- 1 + 2*(a/(1-a))
        real_m <- append(real_m, exp_m)
        old_ms <- append(old_ms, m_suggest)

        new_weights <- get_desh_weights(results_ar1,  p, q, r)
        vars_ar1 <- get_slopes_sides(results_ar1, new_weights, plot_title=sprintf("a:%f, p:%f, q:%f", a, p, q)) #automatically plots RHS 
        new_ms <- append(new_ms, vars_ar1$best_rhs)
        
        ps <- append(ps, p)
        qs <- append(qs, q)
        as <- append(as, a)
        rs <- append(rs, r)
       }
    }
  }
  end_time <- Sys.time()
  print(end_time - start_time)
}
[1] 0.5
[1] 15
Time difference of 0.003687143 secs
[1] 10
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density
[1] 10

[1] 10
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 10
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 10
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 10
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 10
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 10
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 10
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

Time difference of 59.63924 secs
[1] 0.7
[1] 27
Time difference of 0.04512501 secs
[1] 136
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 136
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 136
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 136
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 136
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 136
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 136
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 136
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 136
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

Time difference of 59.47326 secs
[1] 0.9
[1] 120
Time difference of 2.51026 secs
[1] 5995
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 5995
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 5995
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 5995
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 5995
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 5995
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 5995
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 5995
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 5995
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

Time difference of 1.037575 mins
[1] 0.99
[1] 500
Time difference of 7.548814 mins
[1] 119805
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 119805
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 119805
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 119805
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 119805
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 119805
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 119805
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 119805
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

[1] 119805
Warning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true densityWarning: sum(weights) != 1  -- will not get true density

Time difference of 11.9294 mins

# Formatting and storing results
results <- data.frame(old_ms, new_ms, ps, qs, as, real_m)
colnames(results) <- c("Old Method", "New method", "p", "q", "a", "Real Time")
write.csv(results, "./results.csv", row.names = FALSE)
results
# plotting the densities for the best p, q for each a
for(a in poss_as) {
  results_a <- results[results$a == a, ]
  results_a$err_vs_real <- (results_a["Real Time"] - results_a["New method"]) / results_a["Real Time"]
  selected_vals <- results_a[which.min(abs(results_a$err_vs_real$`Real Time`)),]
  acf <- dict[as.character(a)]
  acf_df <- data.frame(acf)
  colnames(acf_df) <- c("acf")
  
  stop_search <- which(acf_df$acf <= 0)[1]
  stop_search <- 1.5 * stop_search
  
  if(a < 0) {
    n <- 3
  } else if(stop_search > 500) {
    stop_search <- 500
  } else {
    n <- 10
  }
  
  y <- get_log(acf_df$acf)[1:stop_search]
  x <- c(0:(stop_search-1))
  lim <- length(y)
  
  results_ar1 <- desh_method(x,y, lim, n)
  print("yo")
  print(results_ar1)
  new_weights <- get_desh_weights(results_ar1,  selected_vals$p, selected_vals$q)
  
  vars_ar1 <- get_slopes_sides(results_ar1, new_weights, plot=TRUE)
  
  
  
  
}
[1] 3
Time difference of 0.0002748966 secs
[1] "yo"
[1] 0
Error in density.default(results$slopes, bw = "nrd", adjust = 0.25, kernel = "gaussian",  : 
  argument 'x' must be numeric
LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKVGhpcyBSIG5vdGVib29rIGNvbmR1Y3RzIHZhcmlvdXMgdGVzdHMgdXNpbmcgdGhlIG1ldGhvZCBvdXRsaW5lZCBpbiBUb3dhcmQgYXV0b21hdGVkIGV4dHJhY3Rpb24gYW5kIGNoYXJhY3Rlcml6YXRpb24gb2Ygc2NhbGluZyByZWdpb25zIGluIGR5bmFtaWNhbCBzeXN0ZW1zLgoKYGBge3J9CiMgRGVwZW5kZW5jaWVzCmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeShMYXBsYWNlc0RlbW9uKQpzb3VyY2UoIi4vZ2VuZXJpYy5SIikKc291cmNlKCIuL2Rlc2hNZXRob2QuUiIpCmBgYAoKYGBge3J9CnN0YXRlc192aXNpdGVkIDwtIGMoKQphY2ZzIDwtIGMoKQpgYGAKCmBgYHtyfQojIFNwZWNpZmljIHRvIGF1dG9yZWdyZXNzaW9uIAphID0gMC45NQpiID0gMgpleHBfbSA8LSAxICsgMiooYS8oMS1hKSkKc3RlcHMgPSAxMDAwMAoKZ2V0X25leHRfc3RlcCA8LSBmdW5jdGlvbihwcmV2KSB7CiAgdyA8LSBybm9ybShuPTEsIG1lYW49MCwgc2Q9MSkKICBuZXh0X3ZhbCA8LSBhICogcHJldiArIGIgKiB3CiAgcmV0dXJuKG5leHRfdmFsKQp9CgpnZXRfaW5pdF9zdGF0ZSA8LSBmdW5jdGlvbiggKSB7CiAgaW5pdF9zZCA8LSBiXjIgLyAoMS1hXjIpCiAgaW5pdF9zdGF0ZSA8LSBybm9ybShuPTEsIG1lYW49MCwgc2Q9aW5pdF9zZCkKICByZXR1cm4oaW5pdF9zdGF0ZSkKfQpgYGAKCmBgYHtyfQojIGV4YW1wbGUgbWV0aG9kIHdpdGggb25lIHNpbXVsYXRpb24gdXNpbmcgb25lIHNldCBvZiBwLCBxCiMgYWxzbyBwbG90cyBoaXN0b2dyYW1zIGhlcmUuIApuIDwtIDEwICMgY2hvc2VuCnAgPC0gMyAjIGNob3NlbgpxIDwtIDMgIyBjaG9zZW4gCnN0ZXBzID0gMTAwMDAKYSA8LSAwLjkKciA8LSAwCgpzdGF0ZXNfdmlzaXRlZCA8LSBhc2ltKGdldF9pbml0X3N0YXRlKCksIGdldF9uZXh0X3N0ZXAsIHN0ZXBzKQphY2ZfZGYgPC0gY19lc3Rfc2VyaWVzKHN0YXRlc192aXNpdGVkKQoKIyBsaW1pdGluZyB0aGUgc2VhcmNoIHRvIGEgZmluaXRlIG51bWJlciwgc2luY2UgdGFrZXMgYSBsb25nIHRpbWUgY29tcHV0YXRpb25hbGx5CnN0b3Bfc2VhcmNoIDwtIHdoaWNoKGFjZl9kZiRhY2YgPD0gMClbMV0Kc3RvcF9zZWFyY2ggPC0gMS41ICogc3RvcF9zZWFyY2gKCnkgPC0gZ2V0X2xvZyhhY2ZfZGYkYWNmKVsxOnN0b3Bfc2VhcmNoXQp4IDwtIGMoMDooc3RvcF9zZWFyY2gtMSkpCmxpbSA8LSBsZW5ndGgoeSkKcmVzdWx0c19hcjEgPC0gZGVzaF9tZXRob2QoeCx5LCBsaW0sIG4pCm5ld193ZWlnaHRzIDwtIGdldF9kZXNoX3dlaWdodHMocmVzdWx0c19hcjEsICBwLCBxLCByKQp2YXJzX2FyMSA8LSBnZXRfc2xvcGVzX3NpZGVzKHJlc3VsdHNfYXIxLCBuZXdfd2VpZ2h0cywgVFJVRSkKdmFyc19hcjEKZXhwX20gPC0gMSArIDIqKGEvKDEtYSkpCnNwcmludGYoIlJlc3VsdGluZyBSSFMgOiAlZiwgZXhwZWN0ZWQgdGltZTogJWYiLCB2YXJzX2FyMSRiZXN0X3JocywgZXhwX20pCmBgYAoKVGhpcyBpcyB0aGUgYmlnIHRlc3QtLSBkb2VzIHRoZSBtZXRob2QgZm9yIG1hbnkgZGlmZmVyZW50IGEncywgcCdzLCBhbmQgcSdzLiAKCmBgYHtyfQojIEdlbmVyYXRlcyBzaW11bGF0aW9ucyBmb3IgZWFjaCBwb3NzaWJsZSBhIHdlIHdhbnQgdG8gdGVzdCwgc3RvcmluZyByZXN1bHRzIGluIGEgZGljdGlvbmFyeS4gCiNwb3NzX2FzIDwtIGMoLTAuOSwgLTAuNywgLTAuNSwgMC41LCAwLjcsIDAuOSwgMC45OSkKcG9zc19hcyA8LSBjKDAuNSwgMC43LCAwLjksIDAuOTkpCnN0ZXBzID0gMTAwMDAKZGljdCA8LSBjKCkKZm9yKGEgaW4gcG9zc19hcykgewogIHN0YXRlc192aXNpdGVkIDwtIGFzaW0oZ2V0X2luaXRfc3RhdGUoKSwgZ2V0X25leHRfc3RlcCwgc3RlcHMpCiAgYWNmX2RmIDwtIGNfZXN0X3NlcmllcyhzdGF0ZXNfdmlzaXRlZCkKICBkaWN0W2FzLmNoYXJhY3RlcihhKV0gPC0gYWNmX2RmCn0KYGBgCgpgYGB7cn0KIyBUaGVuIHdlIHRyeSB0aGUgbWV0aG9kIG9uIGVhY2ggc2ltdWxhdGlvbi4KIyBTdGFydCBieSBmaW5kaW5nIHRoZSBhdXRvY29ycmVsYXRpb24gZnVuY3Rpb25zLiAKcG9zc19wcyA8LSBjKDEsIDIsIDMpCnBvc3NfcXMgPC0gYygxLCAyLCAzKQpwb3NzX3JzIDwtIGMoMCkKCm9sZF9tcyA8LSBjKCkKbmV3X21zIDwtIGMoKQpwcyA8LSBjKCkKcXMgPC0gYygpCmFzIDwtIGMoKQpyZWFsX20gPC0gYygpCnJlc3VsdHNfYXIxcyA8LSBsaXN0KCkKCmZvcihhIGluIHBvc3NfYXMpIHsKICBhY2YgPC0gZGljdFthcy5jaGFyYWN0ZXIoYSldCiAgYWNmX2RmIDwtIGRhdGEuZnJhbWUoYWNmKQogIGNvbG5hbWVzKGFjZl9kZikgPC0gYygiYWNmIikKICBzdG9wX3NlYXJjaCA8LSB3aGljaChhY2ZfZGYkYWNmIDw9IDApWzFdCiAgc3RvcF9zZWFyY2ggPC0gMiAqIHN0b3Bfc2VhcmNoCiAgICAgIAogIGlmKGEgPCAwKSB7CiAgICBuIDwtIDMKICB9IGVsc2UgaWYoc3RvcF9zZWFyY2ggPiA1MDApIHsKICAgIHN0b3Bfc2VhcmNoIDwtIDUwMAogIH0gZWxzZSB7CiAgICBuIDwtIDEwCiAgfQogIAogIGlmKHN0b3Bfc2VhcmNoIDwgNSkgewogICAgc3RvcF9zZWFyY2ggPC0gMTAKICB9CiAgeSA8LSBnZXRfbG9nKGFjZl9kZiRhY2YpWzE6c3RvcF9zZWFyY2hdCiAgeCA8LSBjKDA6KHN0b3Bfc2VhcmNoLTEpKQogIGxpbSA8LSBsZW5ndGgoeSkKICByZXN1bHRzX2FyMSA8LSBkZXNoX21ldGhvZCh4LHksIGxpbSwgbikKICByZXN1bHRzX2FyMXMgPC0gYXBwZW5kKHJlc3VsdHNfYXIxcywgcmVzdWx0c19hcjEpCn0KYGBgCmBgYHtyfQojIFRoZW4gdHJ5IHRoZSBtZXRob2Qgb24gdGhlIGdpdmVuIGF1dG9jb3JyZWxhdGlvbiBmdW5jdGlvbnMuCm9sZF9tcyA8LSBjKCkKbmV3X21zIDwtIGMoKQpwcyA8LSBjKCkKcXMgPC0gYygpCmFzIDwtIGMoKQpyZWFsX20gPC0gYygpCnJzIDwtIGMoKQoKZm9yKGEgaW4gcG9zc19hcykgewogIHByaW50KGEpCiAgc3RhcnRfdGltZSA8LSBTeXMudGltZSgpCiAgYyA8LSA0CiAgTV91cHBlcl9ib3VuZCA8LSA2MDAwCiAgdHJpYW5nbGVfZGYgPC0gZmluZF9tKGFjZl9kZiwgYywgTV91cHBlcl9ib3VuZCkKICBtX3N1Z2dlc3QgPC0gdHJpYW5nbGVfZGZbdHJpYW5nbGVfZGYkcG9zc19NcyA+PSB0cmlhbmdsZV9kZiRjdF9pbnRzLF0KICBtX3N1Z2dlc3QgPC0gbV9zdWdnZXN0JHBvc3NfTXNbMV0KICBtX3N1Z2dlc3QgPC0gdF9pbnQobV9zdWdnZXN0LCBhY2ZfZGYpCiAgCiAgYWNmIDwtIGRpY3RbYXMuY2hhcmFjdGVyKGEpXQogIGFjZl9kZiA8LSBkYXRhLmZyYW1lKGFjZikKICBjb2xuYW1lcyhhY2ZfZGYpIDwtIGMoImFjZiIpCiAgc3RvcF9zZWFyY2ggPC0gd2hpY2goYWNmX2RmJGFjZiA8PSAwKVsxXQogIHN0b3Bfc2VhcmNoIDwtIDEuNSAqIHN0b3Bfc2VhcmNoCiAgICAgIAogIGlmKGEgPCAwKSB7CiAgICBuIDwtIDMKICB9IGVsc2UgaWYoc3RvcF9zZWFyY2ggPiA1MDApIHsKICAgIHN0b3Bfc2VhcmNoIDwtIDUwMAogIH0KICAKICBpZihzdG9wX3NlYXJjaCA8IDUgKSB7CiAgICBzdG9wX3NlYXJjaCA8LSAxMAogIH0KICB5IDwtIGdldF9sb2coYWNmX2RmJGFjZilbMTpzdG9wX3NlYXJjaF0KICB4IDwtIGMoMDooc3RvcF9zZWFyY2gtMSkpCiAgbGltIDwtIGxlbmd0aCh5KQogIHJlc3VsdHNfYXIxIDwtIGRlc2hfbWV0aG9kKHgseSwgbGltLCBuKQogIAogIAogIGZvcihwIGluIHBvc3NfcHMpIHsKICAgIGZvcihxIGluIHBvc3NfcXMpIHsKICAgICAgIGZvcihyIGluIHBvc3NfcnMpIHsKICAgICAgICBleHBfbSA8LSAxICsgMiooYS8oMS1hKSkKICAgICAgICByZWFsX20gPC0gYXBwZW5kKHJlYWxfbSwgZXhwX20pCiAgICAgICAgb2xkX21zIDwtIGFwcGVuZChvbGRfbXMsIG1fc3VnZ2VzdCkKCiAgICAgICAgbmV3X3dlaWdodHMgPC0gZ2V0X2Rlc2hfd2VpZ2h0cyhyZXN1bHRzX2FyMSwgIHAsIHEsIHIpCiAgICAgICAgdmFyc19hcjEgPC0gZ2V0X3Nsb3Blc19zaWRlcyhyZXN1bHRzX2FyMSwgbmV3X3dlaWdodHMsIHBsb3RfdGl0bGU9c3ByaW50ZigiYTolZiwgcDolZiwgcTolZiIsIGEsIHAsIHEpKSAjYXV0b21hdGljYWxseSBwbG90cyBSSFMgCiAgICAgICAgbmV3X21zIDwtIGFwcGVuZChuZXdfbXMsIHZhcnNfYXIxJGJlc3RfcmhzKQogICAgICAgIAogICAgICAgIHBzIDwtIGFwcGVuZChwcywgcCkKICAgICAgICBxcyA8LSBhcHBlbmQocXMsIHEpCiAgICAgICAgYXMgPC0gYXBwZW5kKGFzLCBhKQogICAgICAgIHJzIDwtIGFwcGVuZChycywgcikKICAgICAgIH0KICAgIH0KICB9CiAgZW5kX3RpbWUgPC0gU3lzLnRpbWUoKQogIHByaW50KGVuZF90aW1lIC0gc3RhcnRfdGltZSkKfQoKYGBgCmBgYHtyfQojIEZvcm1hdHRpbmcgYW5kIHN0b3JpbmcgcmVzdWx0cwpyZXN1bHRzIDwtIGRhdGEuZnJhbWUob2xkX21zLCBuZXdfbXMsIHBzLCBxcywgYXMsIHJlYWxfbSkKY29sbmFtZXMocmVzdWx0cykgPC0gYygiT2xkIE1ldGhvZCIsICJOZXcgbWV0aG9kIiwgInAiLCAicSIsICJhIiwgIlJlYWwgVGltZSIpCndyaXRlLmNzdihyZXN1bHRzLCAiLi9yZXN1bHRzLmNzdiIsIHJvdy5uYW1lcyA9IEZBTFNFKQpyZXN1bHRzCmBgYApgYGB7cn0KIyBwbG90dGluZyB0aGUgZGVuc2l0aWVzIGZvciB0aGUgYmVzdCBwLCBxIGZvciBlYWNoIGEKCmBgYApgYGB7cn0KIyBhdHRlbXB0aW5nIGRvd25zYW1wbGluZwphdXRvY29ycnMgPC0gYWNmX2RmCnNlbGVjdGVkIDwtIHNlcSgxLCBucm93KGF1dG9jb3JycyksIDQpCmF1dG9jb3JycyA8LSBhdXRvY29ycnNbc2VsZWN0ZWQsXQphdXRvY29ycnMgPC0gZGF0YS5mcmFtZShhdXRvY29ycnMpCmNvbG5hbWVzKGF1dG9jb3JycykgPC0gYygiYWNmIikKYSA8LSBhCnByaW50KGEpCnN0b3Bfc2VhcmNoIDwtIHdoaWNoKGF1dG9jb3JycyRhY2YgPD0gMClbMV0KICBzdG9wX3NlYXJjaCA8LSAxLjUgKiBzdG9wX3NlYXJjaAogICAgICAKICBpZihhIDwgMCkgewogICAgbiA8LSAzCiAgfSBlbHNlIGlmKHN0b3Bfc2VhcmNoID4gNTAwKSB7CiAgICBzdG9wX3NlYXJjaCA8LSA1MDAKICB9IGVsc2UgewogICAgbiA8LSAxMAogIH0KICAKICBpZihzdG9wX3NlYXJjaCA8IDUgKSB7CiAgICBzdG9wX3NlYXJjaCA8LSAxMAogIH0KICB5IDwtIGdldF9sb2coYXV0b2NvcnJzJGFjZilbMTpzdG9wX3NlYXJjaF0KICB4IDwtIGMoMDooc3RvcF9zZWFyY2gtMSkpCiAgbGltIDwtIGxlbmd0aCh5KQogIHJlc3VsdHNfYXIxIDwtIGRlc2hfbWV0aG9kKHgseSwgbGltLCBuKQpwcmludChyZXN1bHRzX2FyMSkKcDwtMwpxPC0zCm5ld193ZWlnaHRzIDwtIGdldF9kZXNoX3dlaWdodHMocmVzdWx0c19hcjEsICBwLCBxLCByKQp2YXJzX2FyMSA8LSBnZXRfc2xvcGVzX3NpZGVzKHJlc3VsdHNfYXIxLCBuZXdfd2VpZ2h0cywgRkFMU0UpCnZhcnNfYXIxCmBgYA==